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'Jj '. Abstract 

In this paper a general multi-phase-field model is presented which is an ex- 
tension and modification of the model proposed by Folch and Plapp for three 
phase fields [R. Folch and M. Plapp, Phys. Rev. E 72 011602 (2005)] to the 
. arbitrary number of phases. In the model a physical constraint requiring that 

r^ \ the sum of all phase fields in the system is equal to one is resolved by the 

I \ method of Lagrange multipliers. Namely, the thermodynamic driving force is 

^ ■ reduced to its projection on the plane of the constraint. The general model 

O ' functions in a A^- dimensional phase field space were derived which justify 

^ the requirements for the stability of the total free energy functional on dual 

interfaces and hence the absence of "ghost" phases. Furthermore, the case 
^ \ of the different interface energies and mobility parameters on the individual 

0|S ■ interfaces is resolved in a comprehensive manner. It is shown that the static 

equilibrium for three or four phases fulfils Young's law for contact angles 
^ . with high accuracy. Also the model is verified by the quantitative simulation 

■^ ! of the solidification in an Al-Cu-Ni alloy in the case of the four-phase trans- 

formation reaction. We found the way to control the character of new phase 
nucleation using additional terms in free energy functional. 
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1. Introduction 

Multi-phase-field approaches developed in the recent years were practi- 
cally applied to the simulation of the three-phase transformation in eutectic 
and peritectic alloys. The wide range of such realistic microstructure studies 
was reported (see the review article |2| and the references therein). However, 
the investigation of four- and multi-phase transformation reactions is not fully 
covered. The one of reason is the complexity of the models and the emerging 
challenges in terms of the understanding of multi-phase interactions. 

In the last decade the two main concepts of the multi-phase modelling 
have been established. The first one is the multi-phase concept of Steinbach, 
which considers the change of any phase i as the sum of other phase contri- 
butions due to the interaction between the phase i and other phases on all 
individual interfaces, whereas the interaction between more than two phases 
are neglected. Therefore the kinetic of individual interfaces can be considered 
separately with different interface energies and mobility parameters |3|. The 
physical constraint on phase fields, that the sum of all phase fields is equal 
to one is kept in this model automatically. The main works in this area are 
111; Isl, 14], where the main focus is the coupling of the kinetic equations to the 
diffusion equations in the multi-component systems. 

The second concept supposes that the physical constraint on phase fields 
can be resolved by the formal application of the method of Lagrange multi- 
pliers to the free-energy functional which means the geometric projection of 
the driving force vector onto the plane of constraint. The main representa- 
tives of this concept are Folch and Plapp |6| , which developed a qualitative 
phase-field model for three phases with a smooth designed model function 
in the free-energy functional that should ensure the stability of the solution 
and the absence of "ghost" phase on the interfaces between two phases. In 
this model the equations of motion for each of the two solid-liquid interfaces 
can be mapped to the standard phase-field model of single-phase solidifica- 



tion, where the thin-interface asymptotic analysis was applied [15[. Moreover 



it exploits the idea of second order ex pan sion of the free energy functional 



around the equilibrium concentrations |l7l. Il6|). that simplifies the structure 
of the phase-field equations. Later, in the work of Kundin and Siquieri 10 
the total free-energy functional proposed by Folch and Plapp was extended 
by using different thermodynamic factors of phases and was coupled with 
the multi-phase concept of Steinbach. Then, the model was accelerated and 
refined by the inclusion of cross terms in the thermodynamic factor matrix 



in multi-component systems [llj. The same authors [12| apphed the original 
model of Folch and Plapp to investigate the kinetics and morphology of the 
eutectic growth, where the difference in the thermodynamic factors was taken 
into account too. They also showed that the nucleation of new lamellae on 
the boundary of the partner solid phase is rather physical phenomena than 
"artifact" of the model and the probability of this nucleation is defined by 
the the undercooling and the surface tension that is in agreement with the 
nucleation theory. 

The second concept was further exploited in the field of multicomponent 
alloy solidification J8|, |9|. The authors used the Lagrange multiplier and 
wrote the model in common form leaving out of consideration the stability 
requirements. To prevent the existence of the third phase on the individual 
interfaces the authors added the third order term to the obstacle potential 
which hardly can be applied in more general case. Note that, no thin-interface 
analysis of such models is available. 

Using the Steinbach concept and the Lagrange multiplier method an in- 
termediate approach between the first and second concept was proposed, 
which combines the different interface kinetic with the formal account for 



the change of the all phases in the multiple junctions [13|. In the study [14 
it was shown that the reformulated model predicts the angles between the 
phases close to the analytical values except of small derivation in the 3D sim- 
ulations. But by that approach the authors cannot overcome the criticism of 
Folch and Plapp about the instability of the solution for the chosen model 
functions [61]. 

In this paper we use the method of Lagrange multipliers and the idea of 
flatness and stability for model functions suggested by Folch and Plapp |6| 
and then extend it to the iV- dimensional phase- field model (see Section [2]). 
We construct the phase-field model functions step by step, and show the 
implementation in this general model different interface energies and mobil- 
ity parameters on the individual interfaces. The numerical tests presented in 
Section [3] was carried out to show how the model fulfills Young's law, and how 
the mobility parameters influence the dynamic evolution of the correspond- 
ing individual interfaces. In Section H] we verify the model by the qualitative 
simulations of the solidification in an Al-Cu-Ni alloy with the four-phase 
transformation reaction. We first explain the physics of four-phase solidifica- 
tion and write full system of the model equations. We present the numerical 
simulations of the microstructure evolution and show the nucleation effects 
which can be observed. Then we conclude with a summary of our main 



results. 

2. Formulation of general phase-field model 

2.1. Evolution equation for phase fields 

We assume that our physical system can be fully described by A^ phase 
fields Pi G [0, 1], i = 1 . . . N. We identify a vector of phase fields as p = 
{pi, . . . ,Pn), where every phase field Pi means the volume fraction of i-th 
phase. Therefore we demand that in every time our system should follow the 
physical constraint, that is the sum of all phase fields should be equal to one 

N 



Y.P'^ = i- (1) 



i=l 



The total free energy functional of the system is written as 



F= / fdV, (2) 

Jv 

where a total free energy density is expanded into the following terms 

/(p, Vp, c, T) = Kf,{Vp) + HMp) + Up, c, T). (3) 

Here fg sets a free energy cost depending on gradients of phase fields, forcing 
interfaces to have finite width. A constant K has dimension of energy per 
unit length and a constant H has dimension of energy per volume. /;, is a 
dimensionless barrier function, which is analog of double well potential in 
the theory for two phases, fc has dimension of energy per volume and is 
the chemical part of the free energy which depends on concentration vector 
c = (c"^, c^ . . .) and the temperature T. 

For non-conserved phase fields we can write an equation of motion using 
the Ginzburg-Landau equation with Lagrange multiplier. 



, .dpi 1 5F 



1 f6F 1 ^^6F\ , ^ ^^ ,^, 



Here r(p) is a system relaxation time, depending on the phase fields. 

Due to Lagrange multiplier method we produce the projection of driving 
force onto the Gibbs simplex S = {"^^Pi = Ipi & [0, 1]}. That automatically 
ensures ^^ dpt/dt = 0. 



2.2. Construction of model functions with stability and flatness requirement 
The first formulation of stability and flatness requirement for free energy 
model functions was made by Folcli and Plapp in p] as 

-^ L , ni>OVA;, (5) 

and 

SF , 

^=0\fk, (6) 

respectively. In the following we will use requirements which are equivalent 
to Folch and Plapp p, |7| in three dimensional phase-field space but weaker 
in general A^- dimensional space. Namely we require on all interfaces lij = 
{pj = I — Pi, Pkj^i^j = 0} the stability condition 

6^F , 

—^ ^^ >0\/i,j,k, (7) 

and the flatness condition 

and analogical conditions on vertexes Vi = {pi = l,Pk^i = 0} 

6^F , 

— r ^,,^0\/i,k (9) 

6F , 

—- ^, =OVz,A;. (10) 

We start the construction of model functions with choosing the free energy 
gradient term in the form 

/.(vp) = ^E(vp^)'- (11) 

i 

It is easy to check that fg satisfies the flatness conditions (15]) . (ITU]) . 

Then we construct barrier function fh in such a way that we can get 
arbitrary positive interface energy aij for any individual interface lij which 
satisfy our stability and flatness condition (ITl fTO]) . That is 



fb,ij = 2(%U=P. + %U=P,), % = 0'(l-0)'-3(l-0)03 Yl Pfc + 20' Yl Pk ■ 

k^ij kyti,j 

(12) 



It is the polynomial of mininiuni power which follows stability and flatness 
conditions and fb,ijip) = Vp G /fej 7^ lij. Also if p G lij and pi = ip, 
Pj = 1 — if, then fb^ij = (p^{l — (p)^. Therefore we can represent 

fb = ^%/6,ii) (13) 

i<j 

where qij will provide us different surface energies ctjj for every individual 
interface I^j. 

To determine constants qij we should consider the evolution equation of 
phase fields (jl]) on an interface lij. Using fl71 fT0|) it can be shown that only 
i-th and j'-th component of driving force are non zero. Then taking into 
account that all terms fb^ki = on I^j except fb^ij we will see that only /f, j^ 
has contribution in (4). Analogically we can analyze fg term. Therefore 
without chemical free energy term /^ we can write for the static solution of 
eq. © 

dip K d'^ip d(p^(l — (p)^ ^ ,^ ,. 

- %-^ — ^-^^0,t^cx)^ (14) 



dt H dx"^ dip 

d'^ip qijH d(p^{l - (p)^ 



(15) 



dx"^ K dip 

and in equilibrium the total free energy will turn to 

where we use one-dimensional for the sake of simplicity. The solution of 
eq. flTSl) can be expressed as 



x{ip) = I = d(j) for X G (—00, 00), p> G [0, 1]. (17) 



'1/2 y/2Hq,,(j)-\l - (j))^ 
Then we have 

dx \fK 



dip 
dx 


V^v/2^'(l- 


~vY 



(18) 



dip ^2iJg,,v^3(i _ ^Y 

and can express the surface energy as 

/"°° K [^ dip [^ ■ dx 

a., = j_Jdx = -j^ -dip + Hq^,j^ ip\l-ipf-dip 

= ^HKqij [ ^/2ip^{l - ipf dip = ai^HKqij. (19) 
Jo 



Finally qij is defined as 

if^' "^^^^ "^ = T28 
and we can rewrite (fT3l) as 



9*i = U^' where ai = ^^vr. (20) 



/^ = ^]f^E4/Mr (21) 



«<j 



Then constants H and /T can be determined through a model interface width 
W = \/K/H and a maximal interface energy cTmax = maxjj ctjj = aivKH. 
That is K = H^cTmax/'^i and iJ = <Jraa.x/{Wai). Using the above definitions 
we can write barrier function as 



^ 



/^ = E;;^/m.- (22) 



i<j "^max 



V ~ 



Actually in our model we have different numerical interface widths Wi 
l/<Jij. I we change interface energies aij, each particular interface width Wij 
will be changed automatically. We also can change the model interface width 
W in accordance to the need of the numerical method. 

Finally we construct model functions gi{p) which will be used in formu- 
lation of fc- These functions should be equal to one on a vertex Vi and 
on other vertexes. They also should satisfy stability and flatness require- 
ments flTl fTU]) . Then we add a condition for p e lij gi{0, . . .pi, . . . ,pj, . . .0) + 
gi{0, . . .pj, . . . ,pi, . . .0) = 1, which comes from the thin-interface analyses. 

Following the procedure described in Appendix we have found the model 
functions as polynomials of minimal power 



9^ip) = \pI (i5 - 25p, + 15p2 - 3p3 _ 15(1 - p,) Y^p]) ■ 



(23) 



On all individual interfaces lij the functions gi reduces to pf{10 — 15pi + 6pj) 
by taking into account the constraint YlkPk — 1- 

2.3. Evaluation of the smooth function for the mobility parameter 

Here we propose a system mobility parameter T~^{p) which takes constant 
values Tj"^ on any individual interface lij and smoothly varies on the Gibbs 




Figure 1: The illustration plot of the mobility parameterr ^{pi,p2,P3) on the plane pi + 
P2 + P3 = 1, where Tj^^ = 0.5, T23^ = 1, t^^^ = 2. 



simplex S and in neighborhood. Notice that the use of the mobihty parameter 

r~^ allows to consider any immobile interface Jjj with a mobility parameter 

equal to zero t^^ = instead of a relaxation time going to infinity Tij — )• oo. 

Let us to identify a distance Sij between a point inside the Gibbs simplex 



S and an interface lij as 



E 



vl 



(Pi+Pi-l) 



(24) 



Then we can write the mobility parameter of the system as a function of Sjj 
in the form 



T [p) 



i<j 



J ^ij I Z-^^iJ ■ 



(25) 



i<j 



The result is plotted in Fig. 1 where it can be seen that the function (l25l) has 
minimum lines that make it stable if the phase fields changes in the numerical 
simulation outside the Gibbs simplex S. 

The alternative simpler formula which provide the similar simulation re- 
sult can be the following 



-'(p)=E-i>?!';vE 



2 2 
PiPj 



(26) 



i<j 



i<j 



2.4- Evolution equation for concentration fields coupling with equation for 
phase fields 
To evaluate the equations for the concentration fields we consider a multi- 
component system, which contains n chemical components {A, B, . . . ) ex- 
cluding the solvent. From the phase diagram or from the chemical free energy 
functions of phases the equilibrium composition in a phase i can be defined. 
We identify a composition vector Cj and an equilibrium composition vector as 
Ai with components of and Af{T) respectively, a vector of the equilibrium 
chemical free energies of phases with components Bf{Ai,T) = fc^i{Ai,T) 
and a thermodynamic factor matrix of phases Xi with the components 



Note that here we simplified the model presented in [12| by taking the equi- 
librium concentration in the point corresponding to the minimum of free 
energy functions of phases. 

For the definition of the driving force of the phase transformation the 
mixture chemical free energy of a multi-component system can be written 
using the second order Taylor expansion around equilibrium composition 

" yAb 
fc = /- + E V (^^ - ^^'") (^"^ - ^'''") ' (28) 

A,B ^ 

where f^'^ = Yli ^i9i is the mixture equilibrium energy, Bi are chemical 
free energies of phases at the equilibrium composition, c^ = Yli '^^9i ^^^ 
components of a mixture composition, c^'^^ = Y^. Afgi are components of a 
mixture equilibrium composition. 

Then in eq. fl28|) X"^^ are components of a mixture equilibrium thermo- 
dynamic factor matrix, which is defined by the thermodynamic factor matrix 
of phases Xi as 

AT 

X-i = ^Xri^,. (29) 



The derivation of this relation fl29p is given in [12 . 

In the following we will use a mixture diffusion potential vector defined 
in the form of the second order Taylor expansion 

n 
;,^ = 2^X^^.^'?(C^-C^'^''). (30) 



The chemical free energy (I28p gives the thermodynamic part of driving 
force between a phase i and all other phases in the form 



T[p) 



dpi 
dt 



^'(vv^-^Evv) 



dftip) 
dpi 



1 ^ 

-y 



N^ dpk 



+ 



1 ^ 



'iLJ^^'^' 



B. 



(31) 



The diffusion equations for all components have the following form 






V 



J]M^^(p)V/i^-J,1(p) 



(32) 



where M"^^ are the components of the mobility matrix M = D ■ X^^. The 
component of the diffusion matrix are defined as D^^ij)) = Y2i Df^gi with 
Df^ are the terms of the diffusion matrix in a phase i. The values J^ are the 
anti-trapping currents for all components. Then eq. ( !32|) can be transformed 
by the multiplication with X^^ and summation over all components to the 
equation in terms of the chemical potential 



9^ 



^x^^v 



5^M^^V/.^-J,^(p) 



c 



±X--±(^fAf). (33) 



2. 5. Evaluation of the derivatives for model functions Qi 

The full derivatives of the model functions gi according to H] are the fol- 
lowing 






dp, 



dpi N ^—^ dpj 






dp, 



dpi N ^ dpk 



(34) 
(35) 



10 



[(3p.-2)^p5-(1 



where 

I; = ypJ(3|>.-2)}^p^-(l-p,)^(p.-2)l, (36) 

|i = -15p^il -p.)p,, (37) 

5^^ = -l5p^{l-p,)J2p, = -15p^{l-p.)'. (38) 

That yields 

After substitution of (1361) and rearranging we have 

15(2iV-3) 2,, ,2 






R'(i-p.r, (41) 



IEPfc=l 



Pi ( (3Pi - 2) 5^ p2 ^ (3p, + 2) (1 - p2) j (42) 



ap.---^ 2iV ^ ^^^. 

+ 15-pp-p,y-15pp-p,)p,, 

Using eqs. (14T]) and ( 143|) it can be derive that for any individual interface lij 
these derivatives are independent of TV and will be reduced to 15pf{l — pf) 
and — 15p^(l — p]) respectively. 

2.6. Relation to the previous model 

Let us notice that if two polynomials are equal to each other at the Gibbs 
simplex S, then they are equivalent in our model, because they have the same 
derivatives projected at S. Keeping it in mind we found that our functions 
Qi are equivalent to analogical functions suggested by Folch and Plapp in |6[ 
for a 3-phase system. Moreover the derivatives a^|j]pj.=i completely coincide 

11 



with the derivatives of functions gi in Ref. [6|] . The derivatives g^|j]pfe=i 
can be written in the form 

(43) 
{n-2)pi), 



dpi 


EPfc=i 


= 


1 


1) dp, 


EPfc=i 






+ 


15 


)P?a- 






(A^-i; 


a 3-phase 


system 


it will be reduced to 


dpi 


lEPfc=i 


= 


1% 

2% 


EPfe=i 


15 2.. 



Pj)iPk-Pi), (44) 



and is similar to the derivatives of functions gj in Ref. [6|] too. There is only 
the typo in the eq. (3.25) in the Ref. J6|, where instead of -^ is written ^. 
If one will use this equation with typo, then the wrong unexpected behavior 
of Folch and Plapp model can be observed. It can lead to the additional 
nucleation of the third phase (pi) on the individual interface Ijk- 

3. Numerical tests to check Young's law, influence of mobility pa- 
rameters, and absence of ghost phases 

The phase field model should asymptotically convert towards a sharp in- 
terface model when the thickness of interface goes to zero. The corresponding 
sharp-interface model should fulfill the Young's low. The aim of this section 
is to test how does our model fulfill the Young's low in the case of different 
interface energies and mobility parameters. In all simulations here we ne- 
glected chemical potential /c = to test the influence of interface energies 
and mobility parameters and to test whether "ghost" phases exist or not. 
Here all simulations have done in two dimensional coordinate space with the 
fixed boundary conditions. During the tests we have controlled the position 
of the interfaces with mixture of two or 3 phases using the equation @ and 
the position of the points with mixture of 3 phases using equation @. In 
Fig. 2 the initial states for all tests are shown. 

In test 1 we had the system of 64 x 64 discrete points with initial state 
in Fig. 2(b) and investigated evolution of three phases determined in the 
three-dimensional phase-field space with similar interface energies aij = 1 
and similar mobility parameters with Tjj = 1. Then we have carried out 
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Figure 2: Initial configurations for numerical tests, (a) — 64x64 points, (b) — 64x64 
points, geometrically is equivalent to (a), but consists 4 phase fields, (c) — triple point is 
in the center of square, all angles are equal to 120°, (d) — I phase field bounded by II and 
III phase fields in the upper side, the upper boundary equation is y = 64 — tan(7r/6)(a; — 
64)^/64, also II and III phase fields are divided by line x = 64. 



test 2 for four different pliases determined in tlie four- dimensional phase- 
field space with similar initial state Fig. 2(b) and the same parameters. The 
results of the simulation are shown in Fig. 3, which are similar for both tests. 
We found that with a good accuracy both systems have the same dynamic 
evolution that is consistent with expected physics. In the static equilibrium 
which finally is reached we got in 3-phase junction the angles of 120° with 
a high accuracy, which satisfy Young's law. We also can see clear in these 
figures that numerical interface width is constant and we don't have "ghost" 
phases on all individual interfaces. 

In tests 3 and 4 we have simulated 3 different phases defined in the four- 
dimensional phase-field space in the system of 128 x 128 discrete points. The 
system has different interface energies with ratios ai2 '■ ctis '■ cr23 = 0.5 : 0.75 : 
1. In test 3 we have constant mobility parameter r = 1 and in test 4 we have 
different mobility parameters with T12 : Ti^ : T23 = 2 : 1 : 4/3. The initial 
state was taken 712 = 713 = 723 = 120° as shown in Fig. 2(c). According to 
Young's law the respective angles in the static equilibrium should follow the 
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Figure 3: Test 1 and 2. The time evolution of the system presented in Figs. 2(a,b) with 
constant (Ty = 1, r^ = 1. 3 different phases in test 1 or 4 phase in test 2 (first column), 
mixing of 2 or 3 phases (second column) and the mobility parameter T^^(p) (third column) 
at initial, middle and final states. 
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Figure 4: Test 4. The time evolution of the system presented in Fig. 2(c) with ai2 : ctis : 
o'23 = 0.5 : 0.75 : 1, T12 : T13 : T23 = 2 : 1 : 4/3. 3 different phases, mixing of 2 or 3 phases, 
mixing of 3 phases and the mobihty parameter t^^{p) at initial, middle and final states. 



Figure 5: Time evolution of the mixture point of 3 phases for tests 3 and 4. 



ratios sin 712 : sin 713 : sin 723 = ai2 '■ Ci^ : o"23- From this we can calculate 
the angles as 712 = 151°, 713 = 133.5°, 723 = 75.5°. In Fig. 4 the initial, 
middle and final states are shown for test 4. The initial, and final states are 
similar for both tests 3 and 4. The static angels measured are close to the 
theoretical values. It is interesting that we got different numerical interface 
widths Wij, which are inverse proportional to the corresponding interface 
energies ctjj. The different mobility parameters in test 4 have influence only 
on the dynamics of different interfaces, but do not influence the flnal static 
states as it is illustrated by the dependencies in Fig. 5. We also can see in 
Fig. 4 (last column) that the use of eq. [23] give us with a good accuracy 
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{P} 



«j 



^ on every individual interface /, 



ly 



In tests 5,6 (Figs. 6,7) we checked the influence of mobility parameters 
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Figure 6: Test 6. The time evolution of the system presented in Fig. 2(d) with ai2 : (T13 : 
C23 = 1:1:1, T12 : ris : T23 = 2:1: 4/3. 3 different phases, mixing of 2 or 3 phases and 
the relaxation time t at initial, middle and final states. 



on how quickly an interface can be straightened. Namely we used constant 
interface energies aij = 1 in both tests. For test 5 we used constant Ty = 1 
and for test 6 T12 : Tis : T23 = 2 : 1 : 4/3. In Fig. 7 the dependencies of 
the curvatures on the time are shown and we can see again that interfaces 
with higher mobility parameter will move quicker in desired direction, but 
mobility parameters have no influence on the final static states. Due to 
different mobility parameters in dynamic evolution we lose symmetry and 
got it back only in final static state. 
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Figure 7: Time evolution of the curvatures in tests 5 and 6. 




Figure 8: The liquidus surface in Al-reach corner of Al-Cu-Ni phase diagram. 

4. Numerical test for a four-phase reaction 

4.I. Alloy system and model param,eters 

For the numerical test we have chosen a ternary Al-Cu-Ni alloy in the Al 
reach corner of the phase diagram. A zoomed view of the Al reach corner 
in Fig. 8 shows the liquidus surface, the boundary curves and the regions 
where NiAla, Ni2 AI3 and (Al) phases solidify firstly, which we identify as 
a, /3 and 7 phases respectively. Between the a and /3 phases there is the 
peritectic line p^. Between a and 7, /3 and 7 phases there are eutectic lines. 
The three-phase peritectic reaction (ps) and four-phase reactions (f/5) and 
{Uj) are indicated. 

As an example we consider the solidification of an alloys with the initial 
concentration of the liquid of 5.783 at%Ni - 6.278 at%Cu identified as an 
orange point 1 in the phase diagram in Fig. 1. In this alloy crystals of primary 
/3-phase will begin to precipitate at 775°C as the temperature is lowered. 
If cooling continues, the composition of the liquid will change towards the 
boundary curve p^. When the composition reaches the boundary curve at 
730° C (the point 2 in phase diagram), crystals of a-phase will precipitate 
along with crystals of /9-phase. With further cooling, the liquid will change 
its composition along the boundary curve p^ towards the point f/5 while the 
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liquid reacts with crystals of the /3-phase to produce crystals of the a-phase. 
There is the peritectic reaction p^,. 

The main interest is the four-phase reaction in the point f/5 (630°C) where 
crystals of 7 will begin to precipitate along with the /3 phase while the liquid 
reacts with some of the crystals of the a-phase 

L + NiAl3(«) ^ Ni2Al3(/3) + (Al)(7). (45) 

Four various phases coexist in this point. This type of reaction (which in 
many ways is equivalent to the peritectic point on binary diagrams) is known 
as the tributary reaction point (because it looks like a point where two trib- 
utaries of a river meet). The product phase /3 can precipitate along with 
another product phase 7 on the primary phase a, or /3 may be a potent 
nucleant for 7, too. The microstructure formation of such a ternary alloy 
during the directional solidification is of great interest. 

The material and model parameters considered in the simulations are 
listed in Table [H 

Table 1: Material parameters and phase- field model parameters used in the simulation. 



Parameter 


value used 


r (system time scale) 


1 X 10-6 s 


/q (system length scale) 


1.3x10-® m 


Ax/Iq (grid discretization size) 


1 


At/r (time step) 


0.025 


W/Iq (interface width) 


1.1 


D^^ (diffusion in liquid phase) 


1.2 xlO-9 mVs 


D'£'^ (diffusion in liquid phase) 


0.8 xlO-^mVs 


Ds (diffusion in solid phase) 


O.OIDl 


T[/5(4-phase reaction temperature) 


630 °C 


a (surface energy) 


0.30 J/m2 



From the Gibbs free energy functions of phases we estimated the values 
of equilibrium compositions Af, the equilibrium energies Bi and the thermo- 
dynamic factors Xf^ at the temperature 600°C which is below the reaction 
point f/5. The equilibrium parameters of the alloy system are presented in 
Table 2. 



Table 2: Equilibrium parameters at T=600°. Af^\ Af^ and Bi are the normalized equi- 
librium parameters used in the simulations. 



Phase 


Af\ 


Af\ 


B„ 


Ar, 


Af\ 


B. 


xr, 


xf\ 




mol% 


mol% 


J/mol 








J/mol 


J/mol 


L 


7.87 


8.88 


-77231 











2x 10^ 


2 X 10^ 


a 


24.01 


0.1 


-78036 


0.24 





-0.012 


1.0 X 10^ 


1.0 X 10^ 


/3 


17.93 


21.00 


-79537 


0.17 


0.20 


-0.025 


2.0 X 10^ 


2.0 X 10^ 


7 


0.2 


0.2 


-79326 


-0.15 


-0.15 


-0.02 


2.0 X 10^ 


2.0 X 10^ 



4-2. Simulation results 

The four phase reaction were simulated at the constant temperature 600°. 
Equations (13T]) . (l30l) and (l32l) were solved numerically using the Euler method 
in the cubic 2D simulation box of size 200 Ax. The derivatives of the model 
functions gi were calculated according to eqs. (14T]) and (l43l) with A^ = 4. 
The simulations were started with an initial crystal of the a-phase of radius 
12 Ax. After 120 steps a nuclei of the (3 phase was inserted at a random 
site on the solid-liquid boundary of the parent a-phase and after 150 steps a 
nuclei of the 7 phase was inserted in the triple point of the a, /3 and liquid 
phases. 

Results of the evolution of the microstructure at various time steps are 
shown in Fig. 9 (a-d). Three solid phases grow from an initial multi-phase 
nuclei forming the two-phase boundaries. The growth velocity of the a phase 
is slower than the growth velocity of the /3 and 7-phases due to the higher 
Gibbs free energy, so that with increasing time the product phases overgrow 
the crystal of the parent a phase. The chosen model functions serve the 
stability of the solution and the absence of the third phase on individual 
interfaces. The nucleation of 7 phase occurs in the triple point of phases a, 
P and liquid. No nucleation of the third phase on individual interfaces can be 
observed even for the larger undercooling. The lamellar-like structure forms 
by the overgrowing of one phase over other one. 

The time evolution of the Ni and Cu concentration is shown in Fig. 9 (e- 
1). The composition is initially homogen until the evolution of the phase field 
causes the redistribution of the alloy components between the phases. It can 
be seen that the concentration of Ni and Cu in a-phase is larger near the a/7, 
a/P boundaries and smaller at the a-liquid boundary. This phenomenon can 
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be clear explained by the growth of the a-phase in the direction of the liquid- 
phase. The same inhomogeneities in the composition can be observed in the 
/9-phase by the comparison of the /3-liquid boundary and /9-7 boundaries. 

To proof the ability of the iV-phase model to produce the nucleation of 
the third phase we carried out the simulation with an additional nucleation 
term in the kinetic equation: 

where r G [—0.5; 0.5] is a random value. Due to this term a phase i can 
nucleate and grow on the interface between phases j and k if the energetic 
conditions are favorable. The physical meaning is that the nucleation barrier 
can be overcome if the driving force ij,Aj — Bj is large enough. The results 
of the simulation are shown in Fig. 10. It can be observed that a new thin 
(3 and 7-phases form on the /3-liquid and 7-liquid boundaries respectively. 
After increasing time if the a phase is closed and only /3 and 7-phases grow 
the resulting microstructure is similar to the eutectic lamellae structure. The 
same microstructure evolution were produced by the simulation of the eu- 
tectic reaction by means of the 3-phase model of Folch and Plapp. The 



corresponding examples can be found in the works pJl, [18 . 

The comparison of the time evolution of phase fractions for two simu- 
lated case is presented in Fig. 11. The nucleation terms produce the stable 
and uniform growth of /3 and 7 phases with increasing growth velocity. It 
can be shown in the figure that without the nucleation the lamellae of one 
phase overgrow the lamellae of other phase and proceed to grow along the 
solid/liquid boundary where the concentration values are favorable. So far 
the system should wait for the moment when one phase grow enough to go 
around the partner phase. In the case of the additional nucleation terms new 
thin lamellae nucleate at the the solid/liquid boundaries of the partner phase 
immediately after reaching the favorable concentration conditions. The am- 
plitude of the nucleation can be adjusted in accordance to the experimental 
microstructure. 

Notice that we can insert a nuclei arbitrary on a solid/liquid interface 
and the nucleation will occur if the driving force for the nucleation is large 
enough. In the present model the nucleation occurs at the right place at the 
favorable energetic conditions. Moreover the additional terms in the model 
allow to reduce or increase the nucleation barrier, in particular the barrier 
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can be increased in the triple point and prevent the nucleation of the fourth 
phase in this points as it was shown in numerical tests above. 

5. Conclusion 

In this work we have formulated a general phase-field model in A^-dimensional 
phase-field space. The main point of the model is a spatial constructed 
smooth total free energy functional for arbitrary number of phases which 
satisfies the requirements of the stability and flatness on all individual inter- 
faces and vertexes. For this aim the special model functions and responsible 
for the interface energy barrier and for the chemical driving force of the trans- 
formation are proposed which satisfy these requirements and allow to take 
into account the anisotropy of the interface energies and mobility parameters. 

The ability of the model to follow the Young's law and the dynamics of the 
system evolution is tested by the investigation of the phase-field interactions 
at the interfaces and the multiple junctions. It was found that the nucleation 
of new phases can be controlled by the additional terms in the phase-field 
evolution equations. Furthermore, the model is verified by the quantitative 
simulation of the microstructure evolution in an ternary Al-Cu-Ni alloy in 
the presence of the four-phase peritectic-like transformation. 
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Figure 9: Simulated microstructure and concentration fields evolution during isothermal 
four phase reaction in the 2D box without nuleation term. First column represents the 
microstructure at 400 (a), 1000(b) 1800 (c) and 2400 r (d), blue area represents a phase, 
red and yellow areas represent /3 and 7 phases. The corresponding concentration fields of 
Ni and Cu are shown in the second and third columns. 
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Figure 10: Simulated microstructure and concentration fields evolution during isothermal 
four phase reaction in the 2D box with nuleation term ^. First column represents the 
microstructure at 200 (a), 400(b) 700 (c) and 1000 t (d), blue area represents a phase, 
red and yellow areas represent j3 and 7 phases. The corresponding concentration fields of 
Ni and Cu are shown in the second and third columns. 
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Figure 11: Time evolution of the phase fractions for the case without nucleation term (a) 
and with nucleation term (b). 
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